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EXECUTIVE  SUMMARY 


Obj ective 

Develop  high-resolution,  ray-trace  algorithms  for  the  study  of  impulse 
response  functi  ms  in  underwater  sound  channels. 

RESULTS 

1.  A  code  has  been  written  and  documented  illustrating  the  agreement 
with  analytically  tractable  cases,  representative  plots  of  sound  speed 
profiles  and  ray  paths,  and  tabular  output  of  arrival  times,  angles,  and 
transmission  loss. 

2.  Programs,  collectively  called  IMPULSE,  calculated  impulse  response 
functions  for  a  source  and  receiver,  one  of  which  may  be  on  a  boundary  of  the 
duct,  under  the  assumption  of  geometrical  spreading  loss. 

3.  These  programs  allowed  transverse  variation  of  the  refractive  index, 
and  longitudinal  variation  of  the  width  of  the  duct.  Specular  reflection  is 
applied  at  boundaries,  and  bottom  loss  is  included. 

RECOMMENDATIONS 

1.  IMPULSE  should  be  generalized  to  include  range-varying  sound  speed 
profiles  and  absorption  loss. 

2.  Plotting  options  should  be  added  to  permit  plotting  of  arrival 
angles  and  intensity  versus  time.  Ray-path  plots  should  be  altered  to  achieve 
a  set  number  of  points  regardless  of  step  size. 

3.  Applications  of  the  code  in  areas  relevant  to  current  Navy  needs 
should  be  identified.  IMPULSE  should  be  applied  to  experimental  programs  for 
examining  the  predictability  of  actual  ocean  impulse  response  functions. 
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I.  INTRODUCTION 


Itie  FORTRAN  code  presented  here  calculates  exact  {within  machine  word- 
size  and  storaqe  limits)  ray-theoretic  impulse  response  functions.  The 
impulse  response  consists  of  ray  arrival  angles,  travel  times,  and  geometrical 
ray-spreading  transmission  loss.  Ray  paths  are  calculated  by  numerical  inte¬ 
gration  of  a  second-order  ray  equation  of  motion,  and  then  used  to  determine 
flux  tubes  linking  the  source  to  the  receiver,  one  of  which  may  be  on  a  boun¬ 
dary.  The  medium  comprises  1)  a  flat,  lossless,  specularly-reflecting  upper 
boundary  or  "surface,"  2)  a  range-dependent,  lossy,  specularly-reflecting 
lower  boundary  or  "bottom"  3)  a  depth-dependent,  real  index  of  refraction,  and 
4)  an  1/R  attenuation  to  account  for  cylindrical  spreading. 

A  fourth-order  Runge-Kutta  method  is  used  to  integrate 

d2z/dx2  =  -[1  +  (dz/dx)2]c  1  dc/dz ,  (1  ) 

where  z  (x)  is  the  trajectory  and  c  =  c(z)  is  the  phase  velocity.  Equation 
(1)  may  be  viewed  as  the  Euler-Lagrange  equation  implied  by  Fermat's  princi¬ 
ple,  or  simply  a  consequence  of  Snell's  law.  Piecewise  continuous  cubics  are 
used  to  fit  the  index  of  refraction  and  the  bottom  profile,  and  then  used  to 
provide  interpolated  values  between  data  points.  The  bottom  reflection  loss 
is  linearly  interpolated  between  given  data  points. 

The  user  may  enter  arbitrary  increments  for  1  )  the  range  step  size  used 
in  the  numerical  integration,  and  2)  the  uniform  angular  separation  between 
adjacent  rays  in  the  fan  which  originates  at  the  source  and  is  traced  out  to 
the  range  of  the  receiver.  The  limit  angles  in  the  fan  are  also  user  selec¬ 
table.  These  features  are  useful  in  testing  for  the  convergence  of  numerical 
results,  and  for  detailed  studies  of  the  fine  structure  of  arrivals  at  the 
receiver  associated  with  a  particular  angular  region  of  rays  leaving  the 
source . 

This  report  presents  the  computational  features  of  IMPULSE  (ray  tracing, 
surface  and  bottom  reflections,  etc.)  and  a  few  test  cases  of  actual  runs.  It 
describes  the  use  of  the  auxiliary  plotting  routine  provided  and  the  nature  of 
pitfalls  in  applying  IMPULSE  to  problems  with  undersampled  environmental  data 
or  overly  stringent  ray  density  requirements. 
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II.  IMPULSE 


A.  WHY  A  NEW  RAY -TRACE  CODE? 

Frequently,  analysis  of  propagation  requires  knowledge  of  the  basic 
high-frequency,  time-dependent,  multipath  arrival  structure  of  signals  for 
arbitrary  source-receiver  geometries  in  a  particular  analytically  intractable 
waveguide.  Such  results  are  often  good  approximations  to  more  exact 
wavelength-dependent  modal  solutions  and  are  adequate  to  treat  physical  situ¬ 
ations  of  interest,  such  as  underwater  sound  propagation  and  tropospheric 
ducting  of  electromagnetic  waves. 

Many  computer  codes  exist  for  these  or  similar  problems.  Assumptions 
made  in  such  codes  may,  however,  limit  the  generality  of  the  media  one  can 
treat,  or  the  location  of  the  source  or  receiver.  Codes  which  allow  general 
media  may  involve  a  restrictive  amount  of  computation  or  environmental  data  so 
that  the  results  require  considerable  sophisti cation  in  their  use  and  inter¬ 
pretation.  Methods  to  reduce  the  number  of  calculations  may  reduce  the 
accuracy  or  precision  of  the  results  because  the  impulse  response  may  be 
evaluated  only  approximately. 

To  study  the  impulse-response-like  arrivals  of  long-range  underwater 
signals  from  shot-like  sources  (reference  1),  we  recently  developed  a  FORTRAN 
code,  called  IMPULSE,  which  treats  a  depth-dependent  sound  speed  and  a  range- 
dependent  bottom  profile.  Although  limited,  as  discussed  later,  ir  the  types 
of  environments  it  can  treat,  IMPULSE  does  avoid  most  of  the  shortcomings 
mentioned  above.  Since  IMPULSE  may  find  diverse  applications,  this  report  has 
been  prepared  to  give  a  complete  description  of  its  features  and  use.  IMPULSE 
is  not  compared  to  other  codes;  instead,  we  establish  the  agreement  between 
its  calculations  and  analytic  solutions  for  certain  tractable  cases.  Also 
included  is  an  example  using  a  representative,  "real-world"  environment. 

B.  GENERAL  DESCRIPTION 

IMPULSE  computes  almost  exact  solutions  (limited  by  machine  storage)  for 
a  large  class  of  waveguides  in  the  context  of  geometrical  ray  theory.  It 
places  minimal  demands  on  users  and  provides  results  in  a  clear  and  simple 
Fashion.  For  users  with  access  to  DISSPLA  software,  plots  (of  rays  and  the 
curve-fits  to  environmental  data  between  points)  are  provided.  For  other 
installations,  users  may  appreciate  the  option  of  "dumping"  the  information 
used  for  plotting  for  later  access  by  plotting  routines  other  than  DISSPLA *s. 
Compilation  costs  and  basic  storage  needs  are  kept  to  a  minimum  by  writing  all 
ray-path  information  onto  external  files  (if  one  wants  to  get  plots)  and  using 
separately  compiled  plotting  routines.  Since  only  the  point-to-point  features 
of  the  propagation  are  relevant,  no  intermediate  ray-path  points  need  to  be 
stored  in  arrays  in  core. 

The  basic  environment  treated  by  IMPULSE  is  as  follows; 

1.  A  flat,  specularly  reflecting  upper  boundary  or  "surface," 


2.  A  range-dependent  lower  boundary,  or  "bottom,"  described  by  up  to  100 
points 

3.  A  depth-dependent  (transversely-varying)  index  of  refraction,  des¬ 
cribed  by  up  to  100  points,  specified  by  the  phase  velocity  of  the 
wavefronts , 

4.  A  source  from  which  up  to  400  rays  may,  with  uniform  but  arbitrary 
angular  separation,  originate  on  either  boundary, 

5.  A  receiver  which  must  not  be  on,  or  too  near,  a  boundary  (as 
discussed  below),  and 

6.  a  bottom-loss  function  consisting  of  the  loss  in  dB/bounce  at  19 
grazing  angles. 

IMPULSE  traces  rays  from  the  source  to  the  receiver's  range  using  these 
data  and  additional  input  to  govern  the  precision  of  calculations.  It  then 
evaluates  the  standard  expression  (reference  2)  in  ray  theory  for  the  trans¬ 
mission  loss  for  each  pair  of  rays  (or  each  flux-tube)  encompassing  the 
receiver.  This  calculation,  plus  travel  time  and  arrival  angle,  constitutes 
the  impulse  response.  Phase  information  (and  superposition  of  paths)  is  not 
carried  out. 

C.  RECIPROCITY  AND  RECEIVERS  NEAR  BOUNDARIES 

It  is  important  to  note  that  the  usual  expression  for  transmission  loss 
in  geometrical  ray  theory  assumes  a  cylindrically  symmetric  medium  with  the 
source  on  the  axis  of  symmetry.  This  point  bears  directly  on  the  guestion  of 
reciprocity.  In  the  simple  case  of  a  range-invariant  medium,  this  symmetry  is 
trivially  satisfied,  but  if  the  bottom  varies  with  range,  this  assumption  of 
symmetry  automatically  precludes  reciprocity.  The  problem  of  reverse  trans¬ 
mission  (receiver-to-source )  involves  a  different  environment,  namely  one  in 
which  the  receiver  is  on  the  axis  of  symmetry.  Therefore,  no  two-dimensional 
code  with  reciprocity  corresponds  exactly  to  a  physical  problem  if  range 
dependence  is  present.  The  relevance  of  reciprocity  to  the  use  of  IMPULSE  is 
simply  the  desire  to  analyze  a  receiver  on  the  bottom  and,  therefore,  invoke 
the  principle  of  reciprocity:  the  multipath  structure  at  the  receiver  from  a 
suspended  source  is  the  same  as  that  which  would  be  seen  if  the  suspended 
source  were  the  receiver,  and  was  used  as  a  source.  Within  the  tolerance  of 
numerical  errors,  this  type  of  reciprocity  holds  for  IMPULSE  simply  because 
the  ray  tracing  itself  is  carried  out  in  a  two-dimensional  medium  instead  of 
in  a  three-dimensional  medium.  The  cylindrical  spreading  factor  1/R  is 
introduced  purely  artif icially.  In  the  two-dimensional  medium,  a  reduced  form 

2 

of  the  full  three-dimensional  Helmholtz  equation  holds,  that  is,  where  V  = 
2  2  2  2 

3  /3x  +3  /3z  .  The  proof  of  reciprocity  (reference  3)  based  on  this  equa¬ 

tion,  is  still  valid  (with  ray  theory  being  the  high-frequency  limit). 

The  scattering  of  rays  out  of  their  initial  plane  of  propagation  still 
cannot  be  addressed  using  IMPULSE.  Nor,  for  the  reasons  mentioned,  and  since 
IMPULSE  uses  the  standard  geometrical  ray  theory  formula  in  analysis  of  range- 
dependent  three-dimensional  problems  are  the  results  strictly  valid.  Hence, 


the  agreement  between  experiment  and  calculations  will  depend,  in  general,  on 
the  direction  of  propagation  and  bottom  contour. 

For  many  applications,  it  is  reasonable  to  expect  any  discrepancy  to  be 
small.  Analysis  of  the  expected  magnitudes  is  beyond  our  present  aims. 
Similar  comments  naturally  pertain  to  other  ray  or  modal  codes  and  it  is  clear 
that  the  inclusion  of  the  realities  of  three-dimensional  propagation  consti¬ 
tutes  an  enormously  more  complex  family  of  problems,  although  one  which  is 
indeed  seeing  progress  (reference  4). 

The  reason  for  allowing  only  the  source  to  be  on  the  boundary  in  writing 
IMPULSE  is  worth  noting  —  that  there  is  no  basic  difficulty  or  special  coding 
needed  in  tracing  rays  from  the  source.  But,  if  the  receiver  is  on  the 
boundary,  then  the  flux  tubes  encompassing  the  receiver  necessarily  hit  this 
boundary.  Therefore,  the  code  would  have  to  search  for  rays  bracketing  the 
receiver  not  in  depth  (which  is  easy),  but  along  a  generally  irregular  bottom 
(which  is  harder).  It  would  also  have  to  project  the  cross  section  of  the 
flux  tube  normal  to  the  eigenray,  using  rays  at  different  depths  and  ranges, 
instead  of  just  at  different  depths.  For  these  reasons  it  was  chosen  to  allow 
only  the  source  to  be  on  a  boundary  and  to  appeal  to  reciprocity  for  the  cases 
where  this  direction  of  propagation  is  opposite  to  that  in  the  physical 
problem  under  study. 

The  problem  of  a  receiver  too  near  a  boundary  is,  of  course,  still 
present.  Specifically,  suppose  we  have  a  receiver  near  a  boundary  and  trace  a 
fan  of  rays  at  increments  of  2°.  Now  it  may  turn  out  that  no  eigenrays  exist, 
according  to  IMPULSE,  because  no  rays  were  found  to  bracket  the  receiver's 
depth.  It  may  be  clear  by  inspection  that  eigenrays  do  exist.  One  solution 
is  simply  to  determine  roughly  what  rays  at  the  source  are  near  to  eigenrays 
(from  graphs  or  printouts)  and  perform  higher  resolution  ray  tracing  of,  say 
0.2°  increments  in  these  angular  regions.  Since  the  rays  are  now  closer 
together,  it  is  more  likely  that  a  pair  of  them  will  encompass  the  receiver. 
If  not,  the  procedure  can  be  repeated  indefinitely.  Two  pitfalls  can  occur: 
a)  when  the  rays  are  so  close  that  machine  word  size  leads  to  numerical  errors 
in  calculating  cross  sections,  and  ’  )  when  the  receiver  is  exactly  on  the 
boundary  so  that  no  rays  can  possibly  bracket  it  in  depth.  It  is  not  antici¬ 
pated  that  these  cases  are  likely  to  arise  often.  In  such  cases,  it  may  be 
reasonable  to  simply  model  the  problem  as  one  where  the  receiver  is,  in  fact, 
slightly  away  from  the  boundary. 

D.  THE  IMPULSE  ALGORITHM 

This  section  provides  details  on  the  computational  methods  in  IMPULSE. 
For  more  complete  information,  the  user  should  refer  to  the  listing  of  IMPULSE 
in  Appendix  A.  This  listing  may  help  the  reader  familiar  with  the  general 
nature  of  IMPULSE  in  determining  whether  to  carry  out  modifications  in  the 
code  for  more  rays,  refined  bottom  depth  or  bottom  loss,  surface  reflection 
loss,  volume  attenuation,  coherent  summation,  additional  diagnostics,  more 
calculations  (e.g.,  phase  integrals),  or  other  specific  needs  (e.g.,  earth's 
curvature  correction). 

It  may  not  be  essential  to  read  this  section  for  successful  application 
of  IMPULSE.  To  make  best  use  of  the  code,  to  appreciate  problems  that  can 
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arise  in  applications ,  and  to  design  satisfactory  input  data  sets,  it  may  be 
helpful  to  know  some  computational  details.  This  section  also  covers  certain 
points  relating  to  the  accuracy  of  the  code,  introduces  terms  that  might 
otherwise  be  unclear,  and  presents  sign  conventions,  etc.  The  next  section 
describes  the  various  input  variables  and  gives  examples  of  runs. 

It  must  be  stressed  that  IMPULSE  is  basically  a  very  simple  code  with 
highly  restricted,  modest  aims.  The  only  slightly  subtle  aspects  of  IMPULSE, 
perhaps,  are  its  algorithms  for  reflections,  and  its  need  for  adequately 
sampled  data  for  the  index  of  refraction  and  bottom  depth.  The  code  uses 
cubic  polynomials  for  phase  velocity  and  bottom  profiles.  The  piecewise 
continuous  curves  fit  four  data  at  a  time  for  interpolation  between  the  data. 
The  use  of  cubics  leads  to  continuously  varying  functions  with  curvature, 
rather  than  linearly  varying  functions.  (The  use  of  cubics  allows  a  sparser 
sampling  of  data  than  would  linear  segments  to  obtain  equivalent  accuracy  in 
approximation  of  the  original  functions.  Cubics  can  also  lead  to  wild  over¬ 
shooting  between  data  if  the  original  curve  is  undersampled. )  There  are 
several  miscellaneous  details  concerning  the  control  of  precision  of  calcu¬ 
lations,  plotting,  and  diagnostics.  Some  users  will  find  the  use  of  dimen¬ 
sionless  variables  confusing  at  first.  Users  with  electromagnetic  applica¬ 
tions  may  find  the  acoustics  context  disconcerting,  as  we  have  used  the  term 
"sound  speed"  for  phase  velocity  in  the  remainder  of  the  report  as  a  remnant 
of  our  initial  area  of  study.  Additionally,  users  with  applications  other 
than  acoustic  may  need  to  (re)normalize  their  data  in  a  new  system  of  units 
simply  to  conform  to  IMPULSE'S  printing  format  specification,  or,  alterna¬ 
tively,  modify  the  width  of  the  print  fields  in  IMPULSE  to  accommodate  the 
conventional  system  of  units. 

1 .  General  Comments 

IMPULSE  is  written  in  FORTRAN  V.  All  variables  are  double  precision.  It 
consists  of  two  executable  main  routines  and  seven  subroutines.  The  latter 
are  two  routines  for  curve-fitting  the  sound  speed  profile  (SSP)  (or  phase 
velocity)  and  bottom  depth,  a  function  evaluation  (the  right  hand  side  of 
Equation  (1)),  a  cubic  fit  to  the  four  points  on  the  ray  trajectory  if  a 
bottom  reflection  occurs,  a  bottom  reflection  loss  interpolation,  and  two 
routines  for  dunping  the  SSP  and  bottom  depth  curves  (200  points  are  calcu¬ 
lated  and  dumped)  onto  a  file  for  plotting.  The  plot  routine  is  a  separately 
executable  program  from  the  main  routine  The  total  core  storage  required  is 
approximately  1 4K  words  for  the  main  .outine,  and  50K  for  the  plotting 
routine,  after  compilation.  No  external  files  are  required  if  printed  output 
alone  is  desired.  For  plots,  an  external  file  is  needed  with  unit  number  9 
attached  to  it  for  dumping  data.  After  the  run  is  completed  this  file  must  be 
renumbered  as  external  unit  10  for  the  plotting  routine. 

2.  Dimensionless  Variables 

Consistent  length  units  must  be  used  throughout.  If  the  SSP  is  in  m/s, 
depths  must  always  be  in  meters.  It  is  not  permissible  to  mix  feet  and 
meters,  or  meters  and  yards,  etc.,  in  one  run. 
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Input 


The  input  variables  are  described  more  completely  in  the  next  section. 
The  data  consists  of: 

1.  The  SSP  as  (depth,  sound  speed)  pairs. 

2.  A  switch  for  getting  the  SSP  plot. 

3.  The  bottom  depth  data  as  (depth,  range)  pairs. 

4.  A  switch  for  getting  the  bottom  depth  plot. 

5.  The  bottom  loss  in  dB/bounce  at  19  points  in  (bottom  loss,  grazing 
angle)  pairs. 

6.  A  range  step  size. 

7.  The  source  depth. 

8.  The  receiver  depth. 

9.  The  range  to  the  receiver. 

10.  The  angular  spread  of  and  separation  between  rays  leaving  the  source. 

11.  A  tolerance  in  the  precision  of  locating  bottom  reflection  points  in 
range . 

12.  The  number  of  the  ray  whose  trajectory  is  to  be  printed,  along  with 
other  diagnostic  data,  if  desired. 

13.  A  switch  for  getting  plots  of  rays,  and 

14.  The  numbers  of  the  rays  to  be  plotted  if  plotting  is  desired. 

4 .  Ray-Tracing  Method 

Since  IMPULSE  is  basically  short  and  simple,  we  present  its  theory  in 
essentially  complete  form. 

Snell's  law  is  usually  expressed  as  (figure  1) 

n1  sin  41  -  n2  sin  $2.  (2) 

Here,  ^  and  4^  are  the  index  of  refraction  and  the  angle  between  the  ray  and 
the  normal  to  the  surface,  respectively,  in  medium  i. 

If  the  complement  of  4,  is  called  0.  and  n.  defined  as  c  /c,  (where  c  is 

1  11  O  1  O 

a  reference  sound  speed  and  c^  the  local  sound  speed  in  medium  i),  Equation 
(2)  becomes 
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(3) 


C]/cos  0,  =  c2/ccs  6 2  =  cm 

A  coordinate  system  with  dz/dx  =  tan  0  is  now  introduced  and  the  layered 

medium  treated  as  a  continuum;  that  is,  c  =  c(z).  We  define  the  constant 

c(z)/cos  0(z)  as  c  (m  stands  for  maximum),  the  limit-speed  the  ray  can  reach, 
m 

at  which  depth  it  undergoes  total  internal  reflection.  Note  that  z  increases 
in  the  downwards  direction;  0  is  negative  for  an  up-going  ray. 

Squaring  Equation  (3)  and  rearranging,  one  gets 

(dz/dx)2  =  c2  /c2(z)  -  1.  (4) 


Figure  1  The  coordinate  system 


Differentiate  Equation  (4),  obtaining 

2(d2z/dx2) (dz/dx)  =  -2(c2m/c3(z) 3 (dc/dz) (dz/dx) 
or 

d2z/dx2  =  -[c2  /c2(z) ] [1/c(z) ]dc/dz. 
m 

Eliminate  the  first  factor  on  the  right-hand  side  via  Equation  (4)  to  find 

d2z/dx2  =  -(1  +  (dz/dx)2]  c  1  dc/dz,  (5) 

which  can  be  solved  via  Runge-Kutta  methods,  given  c(z). 

Runge-Kutta  numerical  integration  consists,  briefly,  of  computing 
solutions  to  equations  such  as 

d2y/dx2  =  f(y,  dy/dx),  (6) 

given  y(x  ),  at  points  x  +  h,  x  +  2h,  ...  ,  by  evaluating  f  using  weighted 
o  o  o 

sums  for  the  arguments  based  on  the  current  values  of  x,  y,  and  dy/dx.  By  a 
judicious  choice  of  the  weights,  one  can  obtain  accuracy  to  a  high  order  of  h. 
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Reducing  h  shows  whether  or  not  the  numerical  solution  has  converged  on  the 
analytic  solution. 


Specifically,  Equation  (6)  can  be  solved  nunerically  as  follows  (refer¬ 
ence  5):  Let  y'  =  dy/dx,  and  let  y  be  the  numerical  solution  at  x  *  x  + 

n  .  n  o 

nh.  Let  the  numerical  solution  at  xn+1  he  given  by 


y„*1  "  yn  +  h  ,y'n  *  <k1  *  k2  +  k3)/61'  (7al 

y'„„  '  y'n  *  <*.  *  *2  *  *3  *  k4,y6'  ,7b> 

where 

k.  =  hf  (y  ,  y'  ),  (8a) 

i  n  n 

k2  =  hf  (yn  +  1/2  hy'n  +  h^/8,  y'n  +  k,/2),  (8b) 

k_  =  hf  (y  +  1/2  hy*  +  hk  /8,  y*  +  k./2),  (8c) 

J  n  n  i  n  & 

k4  =  hf  (yn  +  hy<n  +  hk3/2,  y'n  +  k^) .  (8d) 


The  difference  between  y  .  and  y(x  .,)  is  proportional  to  h"*.  In  other 

n+i  n+i 

words,  the  Thylor  series  expansion  of  yn+1  agrees,  up  to  and  including  terms 
4 

of  order  h  ,  with  the  Taylor  series  for  y(x  ,).  Proof  of  this  statement  is 

n+i 

algebraically  straightforward,  but  lengthy,  and  is,  therefore,  omitted. 


The  motivations  for  using  the  above  approach  are  its  easy  implementation 
and  the  lack  of  singularities,  or  numerically  ill-conditioned  problems,  at  the 
points  of  total  internal  reflection.  Using  linearly  segmented  SSPs,  where 
rays  are  arcs  of  circles,  is  also  feasible  but  always  biases  the  results?  the 
SSP  linear  segments  lie  inside  every  region  of  curvature.  Linearly  segmented 
SSPs  also  lead  to  more  complex  coding  in  evaluating  the  trajectory  segments. 

For  each  ray  in  the  user-defined  fan  leaving  the  source,  IMPULSE  applies 
the  above  method  to  calculate  the  depth  and  angle  of  the  ray  at  the  range  of 
the  receiver.  The  range  step  size,  h,  is  selected  by  the  user.  Reflections 
from  boundaries  require  special  treatment. 


It  is  stressed  that  each  ray  is  traced  completely  from  the  source  to  the 
range  of  the  receiver  before  the  next  ray  is  traced.  The  endpoint  of  each  ray 
is  stored  as  well  as  the  slope  at  the  endpoint  and  overall  travel  time.  These 
stored  quantities  are  needed  to  compute  the  transmission  loss  and  arrival 
times. 


Unlike  sophisticated  differential  equation-solving  codes,  IMPULSE  does 
not  automatically  estimate  the  error  in  its  numerical  solutions  or  subse¬ 
quently  adjust  the  step  size  to  meet  user-defined  criteria.  IMPULSE  needs  the 
user  to  define  a  step  size  and  uses  that  value  throughout  the  computation. 
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Reflections  at  Boundaries 


IMPULSE  uses  different  methods  to  reflect  rays  at  the  surface  (z  =  o)  and 
bottom.  The  method  for  surface  reflection  actually  affects  the  input  data;  it 
requires  that  the  image  of  the  first  three  data  points  on  sound  speed  below 
the  surface  be  provided  as  input.  That  is,  the  first,  second,  and  third  data 
must  equal  the  seventh,  sixth,  and  fifth  point  except  with  negative  depths, 
the  fourth  datum  is  the  sound  speed  at  the  surface. 

After  a  ray's  depth  is  updated  (advanced  in  range)  the  code  asks  if  the 
new  depth  is  negative.  If  it  is,  the  ray  has  passed  through  the  surface  in 
that  range  step.  Since  the  SSP  is  symmetric  about  the  surface,  all  that  is 
required  then  is  to  reverse  the  sign  of  the  ray's  depth  and  slope  at  the 
updated  point  to  reflect  the  ray  from  the  surface.  In  other  words,  the  image 
of  the  physically  reflected  ray  is  actually  traced  first  and  then  brought  back 
to  physical  space. 

After  a  ray's  depth  is  updated,  and  if  it  has  not  crossed  the  surface, 
the  code  asks  if  the  new  depth  exceeds  the  depth  of  the  bottom  at  the  new 
range.  If  it  does,  a  special  segment  of  code  is  called  into  play,  and  the 
intersection  of  the  local  ray  path  with  the  bottom  is  iteratively  determined. 
The  code  does  not  magically  stop  ray  tracing  at  the  bottom  (nor  at  the 
surface)  but  completes  the  update  using  the  SSP  just  as  if  the  boundaries  were 
absent. 

The  iterative  procedure  works  by  retaining  the  ray's  depth  at  four  points 
—  namely,  the  current  and  updated  values,  plus  the  two  previous  depths.  When 
the  code  finds  that  the  ray  has  passed  through  the  bottom,  it  fits  a  cubic  to 
these  four  points  and  then  determines  the  point  of  intersection  of  this  cubic 
and  the  bottom.  The  iteration  consists  of  working  forward  from  the  current 
range  value  (before  the  ray  depth  exceeds  the  bottom)  in  steps  one-tenth  of 
the  range  step,  until  the  depth  of  the  cubic  exceeds  the  bottom.  Ihen  it  cuts 
the  iteration  step  size  in  half,  reversing  its  sign,  and  works  backwards  until 
the  ray  path  is  above  the  bottom,  and  so  on.  The  procedure  converges  to  a 
point  of  intersection,  since  both  cubics  are  continuous,  but  stops  when  the 
iteration  step  size  becomes  less  than  a  user-specified  value. 

After  the  point  of  intersection  is  determined  to  the  desired  accuracy, 
the  ray  path's  slope  and  the  bottom  slope  are  found.  Specular  reflection  is 
applied  and  bottom  loss  is  calculated  and  added  to  the  previous  bottom  loss  in 
decibel  uni ts . 

Finally,  the  difference  between  the  end  of  the  range  step  and  the  point 
of  intersection  is  calculated  (that  is,  the  distance  the  ray  must  go  to 
complete  the  present  range  step) .  A  second  Runge-Kutta  integration  is  applied 
over  this  distance  and  the  code  resumes  normal  update  procedures.  For 
example,  if  the  bottom  were  hit  one-half  of  the  way  through  the  range  step, 
the  second  Runge-Kutta  integration  would  update  the  ray  over  one-half  a  range 
step  from  the  intersection  point  to  the  end  of  the  range  step. 
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Travel  Time 


Once  the  ray's  depth  is  updated,  IMPULSE  updates  the  accumulated  time 
along  its  path  since  leaving  the  source.  The  arc  length  is  just 

[(dz)2  +  (dx)V/2 

where  dz  and  dx  are  the  change  in  depth  and  the  range  step  size  respectively. 
The  speed  locally  is  taken  to  be  the  average  of  the  sound  speed  at  the  current 
depth  and  updated  depth. 

7.  Transmission  Loss 

IMPULSE  traces  each  ray  within  a  fan  of  rays  leaving  the  source  (at 
user-selected  increments  in  angle  and  between  user-selected  limit  angles)  out 
to  the  range  of  the  receiver,  one  ray  at  a  time. 

As  IMPULSE  traces  a  given  ray,  it  saves  (in  array  elements)  only  the 
current  depth  and  slope,  plus  the  previous  two  depths  (as  variables  to  be 
used,  if  needed,  for  calculating  bottom  reflection  points),  and  accumulated 
travel  time.  Once  a  given  ray  is  completed,  IMPULSE  begins  tracing  the  next 
ray  and  continues  until  all  rays  are  traced  to  the  receiver's  range. 

At  this  point,  the  depth  and  slope  of  each  ray  at  the  receiver  range  are 
stored  in  an  array,  as  well  as  the  "launch"  angle  of  each  ray  at  the  source. 
These  numbers  suffice  to  evaluate  approximately  the  usual  expression 
(reference  2)  for  ray-spreading  transmission  loss  at  the  receiver,  provided  at 
least  one  pair  of  rays,  which  were  adjacent  at  the  source,  bracket  the 
receiver's  depth.  This  requirement  is  the  weakest  constraint  for  approxi¬ 
mating  the  derivative  of  depth  with  respect  to  launch  angle  by  a  ratio  of 
differences  in  depths  and  launch  angles.  If  a  pair  of  rays  not  adjacent  at 
the  source  bracket  the  receiver,  there  is  no  assurance  that  the  depth  is  a 
monotonic  function  of  the  launch  angle  over  the  region  of  values  at  hand. 

Therefore,  there  is  no  reason  to  say  the  ratio  of  differences  is  approximately 

the  derivative. 

Transmission  loss  is  defined  as  the  (decibel)  ratio  of  the  intensity  at 

the  receiver,  and  at  a  unit  distance  from  the  source.  The  formula  for 

transmission  loss  (reference  2)  is 

TL  =  -10  log  [R  '(cos  e  /cos  0  )(dz/d0  )  (9) 

S  R  S  K 

where  R  is  the  receiver's  range  (measured  in  units  of  the  reference  distance 

from  the  source  at  which  the  reference  intensity  is  defined),  0  is  the  angle 

s 

of  the  ray  at  the  source,  0  is  the  angle  at  the  receiver,  and  (dz/d0  )c  is 

R  3  R 

the  derivative  of  the  ray's  depth  with  respect  to  its  launch  angle  at  the 
receiver.  Equation  (9)  is  easily  derived  by  considering  two  rays  leaving  the 
source  at  angles  9g  and  8g  +  respectively.  Their  depth  difference  at  the 

receiver  is  (dz/d$  )  d8  .  Projecting  this  depth  difference  onto  the  cross 
8  R  S 

section  of  the  flux  tube  introduces  a  factor  cos  6  -  The  factor  cos  9  arises 

R  s 


from  the  expression  for  the  solid  angle  subtended  at  the  source  by  the  rays. 
Conservation  of  energy  leads  to  Equation  (9);  that  is,  the  product  of  the 
intensity  and  the  flux  tube's  cross  sectional  area  is  a  constant. 

In  approximating  Equation  (9)  numerically,  both  the  receiver  and  source 
angles  are  taken  to  be  the  average  of  the  angles  of  the  two  rays.  The  deri¬ 
vative  is  just  the  absolute  value  of  the  ratio  of  depth  difference  and  the 
increment  in  source  angle. 

After  Equation  (9)  is  approximately  evaluated,  the  accumulated  decibel 
losses  from  bottom  reflection  for  each  ray  are  added.  Since  the  two  rays 
defining  the  flux  tube  do  not  generally  have  identical  losses,  their  losses 
are  power-averaged. 

8.  The  Use  of  Cubics 

IMPULSE  uses  piecewise  continuous  cubics  to  fit  the  SSP  and  bottom  depth 
data.  The  derivative  in  both  functions  is,  therefore,  discontinuous  at  the 
meeting  point  of  two  cubic  segments.  Ibis  is  not  an  inherent  shortcoming  or 
weakness  in  the  code.  The  Runge-Kutta  formula  always  has  a  well-defined 
gradient  to  work  with,  as  does  the  bottom  reflection  algorithm.  If  spline 
fits  were  used,  an  artificial  curvature  would  arise  in  fitting  a  polynomial  to 
regions  of  linear  variation  following  regions  with  curvature,  simply  because 
the  derivative  would  be  forced  to  be  continuous.  In  our  case,  the  cubic 
segments  lead  to  discontinuous  derivatives.  The  value  of  the  discontinuity 
can  be  easily  made  very  small  simply  by  adequately  sampling  the  curve  being 
fit. 

The  cubic  approximations  can  introduce  large  errors  in  cases  where  the 
cubic  is  forced  to  match  a  very  sharp  bend.  Although  the  curve  always  passes 
through  the  given  data,  it  may  differ  greatly  from  what  one's  eye  determines 
to  be  a  reasonable  curve  between  data  points.  For  example:  if  three  SSP  data 
are  widely  spaced  and  colinear,  and  the  fourth  point  is  only  slightly  below 
the  third  but  with  a  much  different  value,  then  the  resulting  cubic  through 
these  four  points  will  oscillate  wildly  whereas  it  should  be  constant  near  the 
first  three  points. 

The  presence  of  such  unreasonable  fluctuations  is  readily  detected  and 
eliminated  by  means  of  plotting  routines  for  the  curve  fit  representations  of 
the  SSP  and  bottom  depth  functions,  and  the  subsequent  introduction  of  addi¬ 
tional  data  points,  if  needed,  which  lead  the  curve  smoothly  around  the  bend. 
The  added  points  effectively  force  the  cubic  fits  closer  to  the  underlying 
curves.  In  practice,  one  can  ensure  reasonable-looking  curves  by  making 
certain  that  any  regions  of  large  curvature  are  specified  by  at  least  four 
data  points. 

The  cubic  fit  routines  effectively  break  the  data  into  layers  with  four 
points  in  each  layer  and  one  point  in  common  where  layers  meet.  The  first 
layer  consists  of  data  points  (1,2,3, 4)  and  the  second  consists  of  data  points 
(4, 5,6,7)  etc.  Beyond  the  last  compete  layer,  linear  segments  are  used  to 
interpolate  between  data.  Beyond  the  last  data,  the  slope  is  defined  to  be 
zero.  Therefore,  the  index  of  refraction  must  be  specified  by  input  data  to, 
at  least,  the  depth  of  the  bottom  at  its  deepest  to  ensure  that  any  unreal- 


istic  extrapolation  is  avoided.  Similarly,  the  bottom  depth  must  be  specified 
by  input  data  extending  to  at  least  the  range  of  the  receiver  to  avoid 
unrealistic  extrapolations. 

9.  Plotting  Options 


IMPULSE  has  an  independently  executable  plotting  routine  that  uses  the 
D1SSPLA  software.  To  obtain  plots,  a  file  must  be  available  for  IMPULSE  to 
write  onto  as  external  unit  number  9,  and  certain  so-called  input  switch 
variables  (because  the  only  values  recognized  are  0  and  1)  must  be  set  to  1. 
Following  the  IMPULSE  run,  the  code  PLOTS  is  run.  It  will  read  all  needed 
data  from  the  previously  used  file  as  external  unit  number  10  and  requires  no 
other  data.  In  other  words,  the  input  to  IMPULSE  also  defines  each  run  of 
PLOT. 


Plots  are  indispensible  aids  in  debugging  runs  as  well  as  in  interpreting 
the  results.  Plots  show  whether  the  cubic  polynomial  segments  are  good 
representations  of  the  intended  curves,  where  eigenrays  lie,  and  where  bottom 
reflections  occur,  etc. 

IMPULSE  allows  users  to  plot  the  piecewise  continuous  cubic  approximation 
of  the  SSP  at  300  points  uniformly  spaced  between  two  arbitrary,  user- 
specified  depths.  This  is  obtained  by  setting  an  input  variable  called  IPLT1 
to  unity  and  entering  two  depths  between  them.  300  SSP  points  are  calculated 
by  IMPULSE  and  dumped  together  with  the  plot  control  switches. 

IMPULSE  allows  users  to  plot  a  300-point  curve  of  the  cubic  polynomial's 
fit  to  the  bottom  depth  as  a  function  of  range.  In  this  case,  the  points 
automatically  extend  from  the  source  range  (zero)  out  to  the  receiver's  range. 
The  switch  variable  is  IPLT2. 

Finally,  IMPULSE  dumps  ray  trajectory  data,  if  desired,  as  the  rays  are 
traced.  If  PLT3  is  entered  as  unity,  then  the  user  must  specify  how  many  rays 
are  to  be  plotted  and  what  number  position  they  are  in  the  fan  (e.g.,  1,  2,  3, 
4,  5,  6,  7,  8,  9,  10,  30,  31,  32,  33,  34,  35  or  2,  4,  6,  8,  10,  12,  14,  16, 
18,  20,  22,  24,  26,  28,  30,  etc.).  If  IPLT3  is  entered  as  unity,  IMPULSE 

automatically  sets  IPLT2  to  unity  also  so  that  the  trajectories  are  always 
shown  in  relation  to  the  bottom. 

10.  Diagnostics 


IMPULSE  has  a  provision  for  printing  the  evolution  of  up  to  1 00  ray  paths 
in  the  fan.  The  output  also  shows  the  ray's  slope  at  each  point  in  range,  the 
bottom  depth,  and  bottom  slope.  At  bottom  reflections,  the  ray  slope,  bottom 
slope,  and  reflected  slope  are  also  dumped,  plus  the  iteration  history. 

In  addition  to  the  user-directed,  ray-tracing  diagnostic  details,  IMPULSE 
has  a  small  set  of  standard  diagnostic  tests  with  associated  messages  designed 
to  detect  improperly  entered  data  sets,  and  pathological  conditions  arising 
during  ray  tracing,  etc.  These  are: 

a)  A  test  on  the  number  of  SSP  and  bottom  profile  data  which  must  lie 
between  4  and  100  inclusive; 
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b)  A  test  on  bottom-reflected  rays  which  detects  intersection  angles 
less  than  zero  (the  ray  is  approachinq  the  bottom  from  below)  or  greater  than 
90°  (the  ray  is  normal  to  the  bottom  and  would  be  reflected  backwards); 

c)  A  test  on  the  number  of  eigenrays  which  prints  a  message  if  no  flux 
tubes  bracket  the  receiver;  and 

d)  A  test  on  the  bottom  loss  evaluation  which  prints  a  message  if  the 
grazing  angle  (complement  of  the  angle  of  incidence)  is  less  than  the  first 
angle  at  which  bottom  loss  is  given  in  the  input  data  or  is  greater  than  the 
maximum  angle. 

1 1 .  Recommended  Steps  in  Using  IMPULSE 

No  guaranteed  algorithm  exists  to  handle  every  possible  set  of  environ¬ 
mental  data.  In  certain  cases  no  eigenrays  may  be  present.  Other  cases  may 
require  very  careful  searches  for  eigenrays,  or  may  involve  great  sensitivity 
to  small  changes  in  the  environmental  data  or  the  source-receiver  geometry. 
In  the  majority  of  cases,  however,  one  can  expect  to  obtain  reasonable  results 
with  modest  effort.  It  is  still  worthwhile,  nevertheless,  to  list  certain 
basic  steps  a  user  typically  must  take  in  calculating  eigenrays  and  their 
associated  transmission  loss.  (For  a  specific  case  study  in  which  these  steps 
are  followed,  the  reader  is  referred  to  Section  IIIC.)  The  general  steps  are 
as  follows: 

a)  Examination  of  environmental  data.  Prepare  a  data  set  consisting  of 
all  environmental  data  but  with  very  few  rays  requested.  Let  IPLT1  =  1  and 
IPLT2  =  1.  Tftis  creates  plots  of  the  SSP  and  the  bottom  profile.  Execute  the 
adjunct  plotting  program  to  produce  the  plots.  Examine  the  plots  to  see  if 
the  piecewise  continuous  cubic  fits  to  both  functions  are  satisfactory.  If 
not,  add  data  points  in  the  regions  where  the  cubic  disagrees  with  the  intu¬ 
itively  obvious  underlying  curve,  forcing  the  fit  into  better  agreement,  and 
examine  the  new  curves.  When  the  agreement  is  satisfactory,  eigenrays  may  be 
souqht . 

b)  Search  for  eigenrays.  Turn  off  the  switch  variables  for  plotting  the 

SSP  and  the  bottom,  and  set  IPLT3  =  1  to  get  plots  of  the  bottom  and  of  the 
ray  paths.  Select  a  wide  range  of  initial  ray  angles,  a  large  increment 
between  rays  at  the  source,  and  a  large  step  size.  From  the  resulting  ray 
paths'  plots,  or  from  the  printout  of  the  rays'  depths  at  the  receiver  range, 

determine  the  initial  angles  where  eigenrays  are  likely  to  lie.  The  angular 

increment  and  step  size  should  be  large  enough  that  the  runs  made  at  this 

stage  are  not  prohibitively  expensive  and  small  enough  to  give  approximately 

correct  results.  Once  a  family  of  potential  eigenrays  are  found,  they  may  be 
computed  to  greater  precision. 

If  the  search  does  not  suggest  any  likely  candidates,  one  may  slightly 
move  the  receiver  or  source  to  attempt  a  more  favorable  geometry.  Typically, 
this  means  moving  the  receiver  away  from  a  boundary  (see  Section  IIB).  If  it 
proves  difficult  to  find  candidates  after  a  detailed  search,  this  suggests 
that  no  eigenrays  exist.  Physically,  the  signal  reaches  the  receiver  via  an 
evanescent  wave  or  diffraction,  or  via  other  than  great  circle  paths.  In  this 


case,  geometric il  ray  theory  cannot  provide  answers  and  one  must  resort  to 
codes  which  include  wavelength  dependence. 


c)  Detailed  eigenray  computation.  For  the  angular  regions  where  eigen- 
rays  are  expected,  run  IMPULSE  with  a  fine  angular  increment  and  a  fine  range 
step  size.  For  each  pair  of  rays  which  are  adjacent  at  the  source  and  bracket 
the  receiver,  examine  the  transmission  loss  and  travel  time  for  decreasing 
angular  increment,  range  step  size  and,  perhaps,  DXBOT.  When  the  results 
cease  to  vary  when  these  parameters  are  further  reduced,  convergence  is 
obtained.  Studies  of  convergence  are  in  Sections  IIIA  and  IIIB. 

If  results  are  unstable,  examine  the  ray  paths  and  environmental  data  for 
the  presence  of  any  possible  sources  of  this  behavior.  These  may  involve,  for 
example,  the  presence  of  a  knee  in  the  SSP  at  the  exact  depth  of  the  receiver, 
or  a  caustic  surface  near  the  receiver.  Another  potential  source  of  error  can 
be  too  fine  a  step  size;  while  it  has  not  been  seen  in  any  studies  to  date, 
there  is  the  possibility  that  accumulation  of  roundoff  error  over  a  great  many 
steps  can  reduce  the  accuracy  of  IMPULSE.  As  a  crude  rule,  we  suggest  that 
approximately  100  steps  be  allowed  per  loop  length. 

E.  INPUT 

All  input  is  read  under  free-field  format.  Here,  we  assume  that  cards 
are  used,  but  in  the  case  studies,  data  are  from  an  external  file  (unit  8)  as 
in  the  listings  in  Appendix  A.  The  reader  may  refer  to  a  FORTRAN  programmer's 
reference  manual  for  the  system  to  be  used  to  obtain  complete  details  on 
free-field  format.  However,  a  short  description  is  provided  here.  The  free- 
field  format  is  essentially  the  simplest  specification.  If  several  data  are 
on  one  card,  they  only  need  to  be  separated  by  commas.  They  do  not  have  to 
occupy  special  fixed  fields  (i.e.,  columns)  on  the  card.  They  can  be  punched 
in  engineering  or  scientific  notation.  Although  the  variables  read  in  are 
double  precision  they  do  not  have  to  be  represented  as  double  precision  on  the 
cards;  free-field  input  performs  mode  conversion.  If  the  data  read  by  one 
read  statement  extends  over  many  cards,  as  usually  happens  for  SSP  data  and 
bathymetry,  the  end  of  a  card  is  a  legal  delimiter  between  data  items  re¬ 
placing  the  comma.  Therefore,  you  must  not  put  a  comma  after  the  last  item  on 
a  card.  This  would  indicate  to  the  program  that  another  item  is  to  follow  on 
that  card  and  would  result  in  an  execution  error.  If  input  from  another 
device  is  desired,  the  user  must  change  the  unit  nvmber  from  5  to  the  appro¬ 
priate  value  in  the  read  staements  in  the  main  routine,  PULSE.  This  was  done 
for  some  of  our  test  cases . 

Before  describing  the  input  data,  we  first  define  several  FORTRAN 
variables  which  are  read  in  from  cards. 

a)  NPTS1  is  the  number  of  points  in  the  input  sound  speed  (or,  more 
generally,  phase  velocity)  profile  including  the  three  values  above  the  sur¬ 
face  which  create  a  mirror  image  of  those  below  the  surface  (see  section  IID. 5 
on  reflections  at  boundaries). 

b)  IPLT1  is  a  switch  or  control  variable  which,  if  set  equal  to  unity, 
causes  a  set  of  300  points  on  the  SSP  to  be  computed  evenly  spaced  between  two 
user-selected  depths,  after  the  curve-fit  routine  has  been  called,  and  dumped 
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along  with  their  depths  onto  external  unit  9.  If  not  equal  to  one,  IPLT1  has 
no  effect. 

c)  ZZ(I)  is  an  array  of  NPTS1  depths  where  the  input  SSP  data  are  given. 
The  first  three  (if  surface  reflections  are  of  interest),  must  be  the  nega¬ 
tives  of  the  seventh,  sixth,  and  fifth,  respectively  (see  Section  I ID. 5).  I 
must  be  between  4  and  200,  inclusive. 

d)  C(I)  is  an  array  of  sound  speed  values  corresponding  to  depths  ZZ(I). 
(See  remark  at  the  end  of  Section  II. D. 8  regarding  the  extrapolation 
algorithm. ) 

e)  Z1  and  Z2  are  the  depths  between  which  the  SSP  is  computed  from  the 
cubic-fit  routines  if  a  plot  of  the  SSP  is  desired  (that  is,  if  IPLT1  =  1). 
An  array  of  300  values  from  Z1  to  Z2  in  depth  is  dumped  onto  external  unit  9 
if  IPLT  =  1. 

f)  NPTS2  is  the  number  of  bottom  profile  points  in  the  input  data. 

g)  IPLT2  is  a  switch  variable  which,  if  set  equal  to  unity,  causes  a  set 
of  300  depths  on  the  bottom  to  be  computed.  This  is  done  after  the  bottom 
curve  fitting  routine  is  called,  between  zero  range  and  the  receiver  range, 
and  dumped  along  with  their  ranges  onto  external  unit  9.  If  not  equal  to  one, 
IPLT2  has  no  effect. 

h)  RB(1)  is  an  array  of  ranges  at  which  the  bottom  depths  are  input. 

i)  ZB ( I )  is  the  array  of  bottom  depths  corresponding  to  ranges  RB(I). 

j)  BOTLOS ( I )  is  an  array  of  19  bottom  loss  values  in  dB/bounce. 

k)  BOTANG ( I )  is  the  array  of  19  grazing  angles  (in  degrees) 

corresponding  to  BOTLOS,  for  example  0,  5,  ...,  85,  90. 

l)  DX  is  the  range  increment  for  the  Runge-Kutta  routine.  DX  should 
divide  R  (see  below)  an  integer  number  of  times  with  a  small  positive 
remainder. 

m)  ZS  is  the  source  depth. 

n)  ZR  is  the  receiver  depth. 

o)  R  is  the  horizontal  distance  from  source  to  receiver. 

p)  THETA  1  is  the  angle  (in  degrees)  of  the  first  ray  leaving  the  source, 
measured  with  respect  to  horizontal. 

q)  D THETA  is  the  increment  (in  degrees)  in  the  angles  at  which  rays 
leave  the  source. 

r)  THETA2  is  the  angle  of  the  last  ray  to  leave  the  source. 

s)  DXBOT  is  a  cutoff  value  for  iterating  on  the  point  of  intersection  of 
a  local  ray  segment  and  the  bottom. 


t)  ITRACE  is  the  lumber  of  a  ray  whose  trajectory  is  to  be  printed  for 
diagnostic  purposes.  If  ITRACE  *  1,  the  ray  at  THETA  1  is  printed,  for 
example.  If  ITRACE  =  0,  no  trajectories  are  printed. 

u)  IPLT3  is  a  switch  variable.  If  set  equal  to  one,  it  causes  a  set  of 
ray  trajectories  to  be  output  to  unit  9  for  subsequent  plotting.  Otherwise  it 
has  no  effect. 

v)  NPATHS  is  the  number  of  ray  paths  to  be  output  to  unit  9  for  subse¬ 
quent  plotting. 

w)  IPATH(I)  is  an  array  of  the  numbers  of  rays  whose  trajectories  are  to 
be  output  to  unit  9  for  subsequent  plotting. 

Note  in  the  above  description  that  certain  variables  are  irrelevant  if 
certain  of  the  switch  variables  are  not  unity.  If  IPLT1  =  0,  then  Z1  and  Z2 
are  not  needed.  If  IPLT3  =  0,  then  NPATHS  and  IPATH  are  not  needed. 

The  order  in  which  these  input  data  must  appear  on  cards  is  now  des¬ 
cribed.  (IMPULSE  can  be  modified  easily  to  read  from  any  other  external  unit 
than  the  card  reader  simply  by  changing  the  external  unit  number  from  S  to  the 
appropriate  unit  number.)  By  card  number,  the  input  is: 

1.  NPTS1,  IPLT1 

2a.  ZZ(1),  C(1),  ...,  ZZ(NPTSI),  C(NPTS1 ) 

2b.  If  IPLT1  *  1,  Z1,  Z2 

3.  NPTS2,  IPLT2 

4.  RB( 1 ) ,  ZB( 1 ) ,  ...,  RB(NPTS2) ,  ZB(NPTS2) 

5.  BOTLOS  ( 1  ) ,  BOTANGO),  ...,  BOTLOS(19),  BOTANG  (19) 

6.  DX,  ZS,  ZR,  R,  THETA 1 ,  DTHETA,  THETA 2,  DXBOT 

7a.  ITRACE,  IPLT3 

7b.  If  IPLT3  -  1,  NPATHS 

7c.  If  IPLT3  «  1,  IPATH (1),  ...,  IPATH (NPATHS ) 


III. 


CASE  STUDIES 


It  is  essential  to  examine  the  accuracy  of  IMPULSE  as  a  function  of  the 
various  input  variables  to  demonstrate  that  it  does  indeed  converge  to  analy¬ 
tical  answers  with  sufficiently  detailed  input  data.  It  is  also  essential  to 
establish  the  typical  sizes  of  the  "error"  for  typical  values  of  the  number  of 
SSP  data,  range  step  size,  and  ray  angular  spacing.  A  detailed  numerical 
analysis  is  not  presented  here  since  we  are  not  as  interested  in  details  of 
the  various  machine-dependent  tradeoffs  in  truncation  and  roundoff  errors,  as 
in  the  basic  behavior  of  the  code  in  comparison  to  analytic  solutions. 

The  following  two  sections  treat,  respectively,  a  Hirsch  SSP  (reference 
6)  whose  rays  are  sinusoids,  and  a  parabolic  bottom  underlying  a 
constant-speed  medium.  By  suitable  source  and  receiver  positioning,  it  is 
possible  to  find  simple  analytic  expressions  for  transmission  loss  and  signal 
speed  (range  divided  by  travel  time)  for  the  eigenrays.  These  quantities  are 
the  essential  descriptors  of  the  impulse  response  function  for  that  geometry. 
Their  analytic  values  are  compared  to  the  numerical  results  for  separate 
variations  of  NPTS1,  DX,  and  DTHETA  which  are,  respectively,  the  number  of  SSP 
data,  the  range  step-size,  and  the  angle  between  successive  rays  leaving  the 
source . 

It  is  also  essential  to  study  at  least  one  example  of  a  typical  real- 
world  SSP,  including  interaction  with  boundaries.  The  third  section  to  follow 
presents  a  study  of  the  response  at  a  bottom-mounted  hydrophone  to  a  near¬ 
surface  source  for  a  representative  SSP  from  the  Colborn-Wright  data  base  of 
North  Pacific  Ocean  profiles  (reference  7),  using  a  hypothetical  bottom 
contour  and  bottom  loss  function. 

A.  CASE  STUDY  OF  TOTALLY  REFFACTED  RAYS  USING  THE  HIRSCH  PROFILE 

The  SSP  (reference  6) 

c ( z )  -  c  /[I  -  (z-z  )2/a2]1/2  (10) 

o  o 

leads  to  ray  paths  which  are  sinusoidal.  This  is  demonstrated  by  writing 
Equation  (5)  as 

2  2  2  -1 

d  z/dx  =  -(1  +  tan  0)  c  dc/dz 

=  (d/dz) (c2v/2c2)  (11) 

where  cv  *=  c/cos  0  is  the  sound  speed  at  the  "vertex,"  or  point  of  greatest 
off-axis  excursion.  From  Equations  (10)  and  (11), 

d2z/dx2  =  -  (cv/cq)2( z-zo>/a2,  (12) 

which  shows  that  the  ray  obeys  the  same  equation  of  motion  as  a  simple 
harmonic  oscillator. 

For  simplicity,  we  place  both  the  source  and  receiver  on  the  sound  chan¬ 
nel  axis  and  select  the  range  to  equal  one  cycle  distance  of  the  trajectory. 


the 


For  the  initial  condition  where  the  ray  is  at  the  axis,  z(o)  *  z  , 
solution  of  Equation  (12)  is 

z-z  =  A  sin  t (c  /c  )x/a] 

O  V  O 

=  A  sin  (x/a  cos  0Q),  (13) 

where  A  is  a  constant  to  be  determined  and  0Q  is  the  angle  of  the  ray  at  x  = 
0.  Hence,  the  period  of  the  trajectory  (or  loop  length)  is 

x  =  2»a  cos  e  ,  (14) 

P  o 

which  decreases  monotonically  with  increasing  0  .  At  x  *  0,  the  slope  is  tan 

o 

0q.  Differentiating  Equation  (13)  and  setting  x  =  0  therefore  yields 

tan  0  =  A/a  cos  0 

o  o 

or 


A  =  a  sin  0  (15) 

o 

which  shows  that  the  amplitude  of  the  path  (the  maximum  off-axis  excursion) 
increases  as  the  sine  of  the  initial  angle  0Q.  The  final  form  of  the  solution 

is,  from  Equation  (13)  and  Equation  (15), 

z  =  z  +  a  sin  0  sin  (x/a  cos  0  )  .  (16) 

o  o  o 

To  evaluate  the  transmission  loss  in  Equation  (9),  it  is  necessary  to 
have  an  expression  for  the  derivative  of  the  depth  of  a  ray  at  a  fixed  range 
with  respect  to  the  launch  angle  at  the  source.  From  Equation  (16)  this  is, 
at  a  range  of  one  loop  length  (x^), 

,  2 
(dz/d0  )  =  2ira  sin  0  /cos  0  , 

o  xp  o  o 

leading  upon  substitution  into  Equation  (9)  to 

TL  =  20  log  (2«a  sin  0  ).  (17) 

o 

Travel  time  or  signal  speed  is  also  of  interest.  By  definition,  signal 
speed  is  the  range  divided  by  the  travel  time  of  a  ray.  Signal  speed  is  a 
function  of  source  and  receiver  positions.  For  the  case  when  both  the  source 
and  receiver  are  on  the  sound  channel  axis,  it  can  be  shown  that  the  signal 
speed  is 

2 

c  =  2  c  cos  0  /(I  +  cos  0  ).  (18) 

o  o  o 

For  the  numerical  case  study,  we  choose  z  =  z  *  2000,  a  =  3000,  c  = 

so  o 

1480.  The  channel  width  parameter,  a,  and  the  axial  speed,  cq,  are  roughly 
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representative  of  actual  deep  ocean  values  in  meters  and  meters/s,  respec¬ 
tively;  the  SSP  is  quite  unrealistic  for  actual  cases.  The  way  IMPULSE  con¬ 
verges  to  the  analytic  results  for  this  case  should  give  considerable  insight 
into  its  behavior  in  more  general  (analytically  intractable)  cases. 

The  range  is  selected  to  be  one  cycle  distance  of  a  ray  leaving  the 
source  at  15°  depression  angle,  so  that  R  =  x^  =  18207.3  (Equation  (14)).  To 

facilitate  the  execution  of  test  case  runs,  a  special  routine  called  TEST  was 
written  (see  Appendix  B)  which  builds  a  Hirsch  SSP  data  file  for  IMPULSE  but 
in  itself  needs  only  general  instructions  such  as  the  number  of  SSP  data,  zq, 
a,  R,  DX,  and  plotting  switches. 

Figures  2  and  3  show  the  SSP  and  a  set  of  ray  paths  for  this  case.  In 
figures  4  through  8,  we  present  a  comparison  of  the  numerical  results  and  the 
analytic  results.  The  analytic  result  for  transmission  loss  (Equation  (17)), 
is  TL  =  73.8  dB  and  for  signal  speed  (Equation  (18)),  it  is  c  •=  1479.1. 

Figure  4  shows  that,  except  for  residual  errors  from  other  parameters, 
the  code  converges  to  the  analytic  value  for  TL  if  30  points  are  input  to 
describe  the  SSP.  The  anomaly  at  NPTS1  «  14  is  not  understood,  but, 
presumably,  it  arises  because  a  layer  in  the  piecewise  cubic  fit  to  the  SSP 
ends  at  the  channel  axis  and  another  begins.  As  noted  in  the  figure  caption, 
DX  -  303.45  and  DTHETA  -  0.1°.  With  this  value  of  DX,  IMPULSE  makeB  60  range 
steps  to  reach  a  range  of  1 8207 ,  which  is  0.3  away  from  the  range  where 
analytic  results  are  evaluated.  Therefore,  we  have  60  range  steps  per  loop 
length. 

Figure  5  shows  that  the  numerical  resultB  for  TL  converge  to  the  analytic 
results  when  as  few  as  8  range  steps  are  taken  per  loop  length.  In  this  case, 
NPTS1  -  100  and  DTHETA  -  0.1*. 

Figure  6  presents  a  study  of  the  sensitivity  of  signal  speed  results  to 
DX.  The  same  parameters  are  used  as  for  figure  5.  Convergence  is  obtained 
for  DX  less  than  approximately  0.025  xr  (i.e.,  40  range  steps  per  cycle). 

Figures  7  and  8  show  the  effect  of  varying  the  spacing  of  rays  at  the 
source,  DTHETA.  The  plots  show  a  pair  of  sudden  jumps  in  the  numerical 
results.  Examination  of  the  output  from  IMPULSE  (not  shown)  in  the  neighbor¬ 
hood  of  these  jumps  explains  the  behavior.  Since  IMPULSE  searches  for  rays 
which  bracket  the  receiver,  it  is  very  sensitive  to  slight  changes  in  inputs 
which  can  switch  the  particular  rays  bracketing  the  receiver.  At  the  position 
of  the  jumps,  the  detailed  output  from  IMPULSE  shows  that  the  pair  of  rays 
found  by  IMPULSE  does  indeed  hop,  so  that  the  TL  or  signal  speed  is  basically 
calculated  on  rays  largely  above  the  receiver  and  then  largely  below  the 
receiver.  In  other  words,  the  receiver  lies  on  the  extreme  edge  of  a  ray 
pair.  The  precise  position  of  such  jumps  depends  on  how  one  defines  the  fan 
of  rays  (the  limit  angles  91  and  62)  and  their  spacing.  Convergence  to  analy¬ 
tic  results  (to  within  residual  errors)  is  obtained  for  angular  spacing  less 
than  or  approximately  equal  to  0.2°. 
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Figure  2.  The  SSP  of  equation  1 0  is  shown  for  the  case 
zQ  =  2000,  a  =  3000 
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Figure  3.  Ray  paths  are  shown  for  01  =  5°,  02  =  20°,  DO  =  1 .0  °, 
DX  =  364. 1 4,  and  NPTS1  =  100. 
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Figure  4.  The  numerical  values  (dots)  are  compared  to  the  analytic  value  (dashed 
line)  of  transmission  loss  for  various  values  of  the  number  of  SSP  data.  The 
anomaly  at  NPTS1  =  14  is  not  understood.  The  step  size  was  DX  =  303.45 
(R/DX  ~  60)  DTHETA  =  0.1  °  Circles  represent  variations  of  DX  for  fixed 
NPTS1  =  14.  Residual  errors  probably  arise  from  DX  and  DTHETA. 
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Figure  5.  The  numerical  transmission  loss  (dots)  is  compared  to  the  analytic 
result  (dashed  line)  for  several  values  of  DX,  the  range  step  size,  corresponding 
to  4,  5,  ...,  10  steps  in  a  cycle  of  the  trajectory.  For  these  results,  NPTS1  =  100 
00  =  0.1° 


Figure  6.  The  numerical  values  of  signal  speed  (range  divided  by  travel  time)  are 
compared  to  the  analytic  result  (dashed  line)  for  varying  range  step  size  in  the 
Runge-Kutta  solution  for  the  ray  path.  For  these  runs  NPTS1  =  100  and 
DTHETA  =  01°. 


Figure  7.  The  numerical  values  of  transmission  loss  (dots)  are  compared  to  the 
analytic  result  (dashed  line)  for  variations  in  DTHETA.  (The  jumps  at  1 0  and  3° 
result  from  the  combination  of  DTHETA  and  the  angular  range  limits  THETA1 
and  THETA2,  which  leads  to  switches  from  eigenray  pairs  largely  above  to  largely 
below  the  receiver.) 
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DTHETA  (deg) 

Figure  8.  Numerical  results  for  signal  speed  (dots)  are 
compared  to  the  analytic  result  (dashed  line)  for  various  values 
of  DTHETA,  the  angular  increment  between  adjacent  rays 
leaving  the  source.  For  these  results,  NPTS1  =100  and 
DX  =  364.14.  (Jumps  at  1  °  and  3°  arise  from  the  choice  of 
DTHETA  and  the  limiting  angles  of  the  fan  of  rays  at 
the  source.) 


It  is  worthwhile  to  point  out  that  these  results  are  specific  to  the 
particular  SSP  studied.  For  other  cases  one  should  generally  anticipate 
different  sensitivities,  even  though  the  required  precision  in  input  may  be 
roughly  similar  to  our  case  study. 

B.  CASE  STUDY  OF  BOTTOM -REFLECTED  RAYS  USING  A  PARABOLIC  BOTTOM 

The  preceding  section  presented  a  study  of  the  accuracy  of  IMPULSE  for 
rays  which  do  not  hit  the  boundaries,  but,  instead,  are  totally  refracted. 
Another  important  question  of  accuracy  concerns  bottom-reflected  paths.  A 
convenient  environment  for  studying  this  class  of  problems  is  one  in  which 
rays  are  straight  line  segments  while  the  bottom  forms  a  parabolic  mirror  with 
the  source  at  the  focus.  In  this  case,  all  rays  hitting  the  bottom  are 
reflected  to  travel  horizontally  (in  the  analytic  solution).  We  develop 
formulas  for  transmission  loss  in  this  case,  compare  numerical  results  with 
analytic  solutions,  and  study  the  angles  of  the  reflected  rays. 

Let  the  bottom  depth  z  ,  as  a  function  of  range,  be 

D 

z  (x)  -  z  +  (a2  +  2ax)1^2  ,  (19) 

D  O 

so  that 

x(z  )  =  [ (z  -z  )2-a2)/2a.  (20) 

B  BO 

Equation  (20)  represents  a  parabola  with  a  focus  at  range  x  =  0  and  depth  z  . 
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To  calculate  the  tranemision  loss  for  bottom-bounce  rays,  it  is  necessary 
to  obtain  a  relation  between  the  infinitesimal  change  in  the  launch  angle  and 
the  infinitesimal  change  in  the  ray's  depth  at  the  receiver's  range.  Figure  9 
shows  a  construction  useful  in  obtaining  the  desired  relation.  Let  X  and  Y  be 
the  points  of  intersection  of  the  rays  launched  at  6fi  and  08  +  40g>  respec¬ 
tively,  with  the  bottom.  Then,  with  0_  being  the  bottom  angle, 

D 

XY  sin  (6e  -  0g)  -  Ir  A0g  ,  (21) 

XY  sin  9  -  Az, 

D 

where  SR  is  the  "slant  range,"  or  distance  from  the  source  to  X.  Since  AO^  is 

very  small,  the  value  of  SR  is  the  same  to  point  Y  as  to  X  within  terms  of 
order  A6g.  We  can  ignore  the  higher  order  correction. 


ray  depth  at  the  receiver’s  range  with  respect  to  the 
launch  angle. 


Then 


SR 


r  +  a 


(22) 


where  r  is  the  range  of  the  point  where  the  ray  hits  the  bottom. 

From  Equations  (21)  and  (22), 

dz/d0_  =  (r  +  a)  sin  6  /sin  (0  -0  ).  (23) 

S  BSB 

To  evaluate  Equation  (23),  values  are  needed  for  r  and  0R.  The  range  is 

found  directly  from  Equation  (19)  for  the  bottom  depth  as  a  function  of  ran'  ^ 
which  can  be  inverted  for  range  as  a  function  of  depth  (Eq.  20)).  0R  is  foui.J 

by  differentiating  the  expression  for  depth  with  respect  to  range.  Since  the 
rays  become  horizontal  after  reflection,  the  bottom  depth  at  the  point  of 
intersection  is  the  same  as  the  receiver  depth,  and,  hence,  at  point  X,  the 
range  is 

x  *  [(zR  -  Zg)2  ~  a2] /2a .  (24) 

Finally  we  have 

dz  /dx  =  tan  0_  -  a/[a2  +  2ax ]^2 

B  3 

or 

sin  0  *=  a/[2a2  +  2ax]1y^2  .  (25) 

B 

To  evaluate  the  transmission  loss,  0„  is  alBO  needed  and  is  obtained  from 

S 

tan  6  ”  (z_  -  z  )/r . 

S  E  S 

For  our  numerical  study  we  choose  z  “  1000,  z_  ■  3000,  R  •  5000,  a  = 

OK 

1000.  An  example  of  the  bottom  profile  and  the  ray  paths  for  these  parameters 
is  given  in  Figure  10.  As  in  the  Hirsch  SPP  study,  we  have  used  a  special 
routine  TESTB  that  creates  a  data  file  for  use  by  IMPULSE  which  requires  only 
a  small  number  of  input  data  such  as  NPTS2  (the  number  of  bottom  profile  data 
to  be  created),  a,  DXBOT ,  zs>  zR,  DX,  R,  01,  D0 ,  02,  and  plotting  switches.  A 

listing  of  TESTB  is  presented  in  Appendix  B. 

For  this  example,  the  range  at  which  the  eigenray  strikes  the  bottom  is 
1500.  The  bottom  slope  at  this  point  is  0.5',  and,  hence,  0_  “  26.57'.  The 

eigenray  leaves  the  source  with  slope  1.333,  and,  hence,  eg  -  53.13*.  The 

transmission  loss  is  TL  «  73.2  dB. 

Figure  11  shows  a  comparison  between  numerical  results  and  the  analytic 
result  for  TL  as  a  function  of  the  number  of  points  used  to  define  the  bottom 
profile . 
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Figure  10.  Ray  paths  and  bottom  depth  are  shown  for  a  parabolic-mirror  bottom. 
For  these  results,  NPT2  =  50.  a  =  1000,  DXBOT  =  1,  ZS  =  1000.  DX  =  50,  =  40°. 

D0  =  2°,  and  62  -  60° 


Figure  12  shows  the  scatter  in  ray  angles  reflected  from  the  bottom  and 
as  a  function  of  NPTS2.  In  both  the  TL  and  angle  calculations,  convergence  is 
relatively  good  for  NPTS2  greater  than  25. 

C.  CASE  STUDY  USING  A  REPRESENTATIVE  SOUND  SPEED  PROFILE  IN  THE  NORTH 

PACIFIC  OCEAN 

This  section  presents  a  typical  example  of  the  use  of  IMPULSE  to  cal¬ 
culate  the  impulse  response  for  a  bottom -mounted  receiver  and  a  near-surface 
source  in  an  SSP  typical  of  the  North  Pacific  Ocean.  Also  included  are  a 
(hypothetical)  bottom  loss  function  and  bottom-depth  function.  The  primary 
aim  here  is  not  to  obtain  extremely  high  accuracy,  but  only  to  demonstrate, 
for  a  specific  numerical  case,  how  one  would  apply  the  general  procedure  for 
using  IMPULSE  as  discussed  in  Section  IID.  This  includes  typical  problems 
often  encountered  in  practice  such  as  obtaining  reasonable  fits  to  the  SSP  and 
locating  eigenrays. 

This  example  calculation  shows  that  an  uncertainty  in  the  end  results  is 
unavoidable  because  of  the  finite  number  of  environmental  data,  and  the 
ensuing  freedom  allowed  in  generating  reasonable  curves  to  interpolate  between 
the  given  data.  A  personal  judgement  must,  ultimately,  be  introduced  in  this 
area,  The  example  also  shows  a  characteristic  of  long-range  undersea  ray 
paths,  namely  great  sensitivity  of  the  endpoint  depth  and  angle  to  the  initial 
angle. 
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Figure  1 1 .  For  the  parabolic  mirror  bottom,  the  transmission  loss  from  IMPULSE 
is  compared  to  the  analytic  result  (dashed  line)  for  different  numbers  of  bottom 
depth  data  NPTS2.  In  these  runs,  DX  =  200,  DTHETA  =  1  °,  R  =  5000,  zs  =  1000, 
zR  =  3000,  DXBOT  =  10,  a  =  1000,  and  the  sound  speed  is  constant. 
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Figure  12.  The  angles  of  computed  rays  at  their  endpoints  after  bottom  reflection 
are  shown  as  a  function  of  NPTS2.  For  these  cases,  rays  left  the  source  at  50°, 
51  °, ...  60°  and  all  other  parameters  were  the  same  as  in  figure  1 1  and  in 
the  text 


The  first  step  is  to  examine  the  environmental  data  and,  if  necessary, 
add  or  delete  data  so  as  to  obtain  reasonable  curve  fits  to  the  data.  The  SSP 
chosen  for  this  case  study  is  from  the  Colborn-Wright  data  base  (references 
7,8),  province  19,  Fall,  in  the  North  Pacific  Ocean.  Figure  13  shows  the 
input  data  for  a  trial  run.  There  are  38  SSP  data  pairs.  The  SSP  plot  data 
is  dumped  (onto  external  unit  number  9)  between  0  and  1000  m.  A  set  of  19 
bottom  depths  is  specified  out  to  range  200  km.  The  bottom  loss  is  purely 
hypothetical.  The  source  depth  is  1000  m,  the  same  as  the  bottom  depth,  a 
zero  range.  The  range  step  size  is  1000  m  and  the  receiver  is  at  a  depth  of 
2000  m  and  range  of  200  km.  Therefore,  in  this  case,  rays  actually  are 
traced  from  receiver  to  source  and  reciprocity  invoked.  A  set  of  six  rays  is 
traced  from  -10*  to  0*  in  2®  steps.  Bottom  reflections  are  found  to  within  5 
m  in  range.  No  printed  trajectories  are  requested.  Plots  are  requested  for 
all  six  rays. 
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Figure  1 3.  The  input  data  are  shown  for  the  first  trial  run 
of  IMPULSE  using  an  SSP  representative  of  the  North  Pacific  Ocean. 


As  a  point  of  interest,  note  that,  before  this  run  was  made,  the 
necessary  file  for  storing  plotting  data  was  created  and  the  external  unit 
number  9  attached  to  it.  The  input  data  was  stored  on  a  file  with  external 
unit  number  8  and  IMPULSE  was  modified  so  as  to  read  from  unit  8  rather  than 
from  5  (the  card  reader). 

Figure  14  shows  the  printed  output  from  IMPULSE  for  the  trial  run. 
Although  an  eigenray  pair  i3  found,  the  calculation  is  not  at  a  very  high 
precision.  It  is  to  be  expected  that  further  runs  with  greater  precision 
(smaller  range  step)  will  move  the  rays  to  new  positions. 
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Figure  15  shows  the  SSP  plot  from  PLOTS,  the  plotting  adjunct  which  uses 
the  output  file  on  unit  9  as  input.  (This  file  is  attached  to  unit  10  before 
running  PLOTS.)  The  curve  fit  is  unsatisfactory  because  of  the  spurious  local 
maximum  at  about  200  m  depth.  There  is  no  unique  procedure  for  generating  a 
better  SSP  data  set  and  one  must  use  personal  judgement  in  adding  or  deleting 
data. 
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Figure  1 5  The  SSP  as  fit  by  the  piece-wise  continuous  cubics  of  IMPULSE.  Note  the 
spurious  local  maximum  near  the  200  m  depth. 

The  solution  of  adding  data  was  employed  in  this  test  case.  On  a  sepa¬ 
rate  plot  (not  shown),  the  initial  SSP  data  were  plotted  to  a  depth  of  500  m. 
Then  a  smooth  curve  was  sketched  between  data.  Where  Figure  15  shows  unrea¬ 
sonable  curvature,  new  data,  lying  on  this  interpolated  curve,  were  inserted. 
The  input  data  file  was  augmented  with  these  data  at  depths  of  40,  175,  225, 
275,  350,  and  450  m.  Figure  16  shows  the  new  data  file  and  figure  17,  the  new 
SSP,  which  is  satisfactory  for  present  purposes. 

A  procedure  exactly  analogous  to  that  for  the  SSP  data  revision  can  be 
applied  to  smooth  any  undesirable  curve-fitting  artifacts  in  the  bottom  depth 
function.  For  this  test  case,  it  turns  out  that  the  curve  fit  is  reasonably 
free  of  undesirable  features. 

The  revised  printout  from  IMPULSE  is  shown  in  Figure  18  —  the  eigenrays 
are  gone. 

The  next  step  is  to  locate  eigenrays.  Figure  19  shows  the  ray  paths 
obtained  from  the  second  run  with  a  revised  SSP. 

As  a  point  of  interest,  the  user  should  be  aware  that  plots  showing  rays 
reflecting  from  the  bottom  often  will  not  show  the  actual  point  of  contact. 
In  other  words,  the  ray  will  appear  to  bounce  from  a  point  above  the  bottom. 
This  does  not  imply  errors  have  occurred  in  the  ray-tracing  or  plotting 
routines.  It  actually  arises  because  the  plots  of  ray  paths  are  made  from 
points  which  are  spaced  by  DX  in  range.  The  plotting  routine,  of  course,  has 
no  means  of  knowing  where  the  ray  went  between  these  points.  Consequently, 
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two  points  before,  and  after,  a  bottom  reflection  will  just  be  connected 
directly  by  the  plotting  software,  without  any  line  segment  touching  the 
bottom.  The  same  artifact  will  also  be  seen  for  surface  reflections. 

The  sound  speed  at  the  source,  and  at  the  receiver,  are  useful  in  sug¬ 
gesting  where  eigenrays  may  lie.  From  Snell's  law,  the  ratio  of  sound  speeds 
can  be  used  to  determine  a  limiting  ray  angle,  at  the  source,  which  must  be 
exceeded  simply  in  order  for  the  ray  to  reach  the  receiver.  In  our  case,  the 
sound  speed  at  the  source  is  1481.3  m/s  and  at  the  receiver,  1492  m/s, 
yielding  a  limiting  ray  angle  of  approximately  7°.  Since  the  bottom  slope  at 
the  source  is  roughly  1  ° ,  any  eigenrays  must  leave  the  source  in  the  upward 
direction.  Therefore,  as  an  initial  fan  of  rays  suggests  eigenray  locations, 
we  take  initial  and  final  angles  of  -30°  and  -5°,  respectively,  and  a  step 
size  of  5°  between  rays.  The  range  step  size  is  chosen  to  be  100  m,  corre¬ 
sponding  to  2000  steps  from  source  to  receiver.  Figure  20  shows  the  ray  paths 
and  Figure  21  the  printout  for  this  case.  Three  eigenray  pairs  or  flux-tubes 
are  identified. 

There  is  no  guarantee  that  the  eigenrays  identified  at  this  stage  are 
valid,  because  the  step  size  has  not  been  reduced  to  the  point  where  further 
reduction  causes  no  change.  In  a  set  of  tests  (not  shown),  we  find  results 
stabilize  for  values  of  DX  around  10  m.  This  is  a  remarkably  small  value, 
corresponding  to  20,000  steps  from  source  to  receiver.  In  general,  each 
environment  must  be  evaluated  in  its  own  context  and  some  cases  will  not 
require  as  many  calculations.  But,  for  bottom  profiles  with  sharp  angles  or 
other  special  cases,  one  can  expect  extreme  sensitivity  of  ray  endpoints  to 
small  changes  in  launch  angle. 

The  remaining  flux  tubes,  after  results  have  converged,  are  clustered 
around  -24.8°  and  -24.4*  at  the  source.  Figure  22  shows  this  in  printout 
form.  For  these  rays,  it  is  not  desirable  to  obtain  plots  since  such  a  huge 
number  of  points  would  be  involved  with  the  existing  routines.  It  would  be 
more  straightforward  to  modify  the  plot  routine  so  that  each  path  would  be 
described  by  200  points  (e.g.,  by  throwing  out  points),  but  this  has  not,  as 
yet,  been  accomplished. 
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Figure  18.  The  revised  IMPULSE  printout.  Note  that  no  eigenrays  are  found. 


Figure  1 9.  Ray  paths  and  bottom  profile  are  shown  for  the  test  case  —  for  launch  angles  from  -10°  to  0° 
in  2°  steps.  Bottom  reflected  rays  may  appear  to  reflect  from  points  above  the  bottom. 
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Figure  21 .  The  printout  for  a  fan  of  rays  from  -30°  to  -5°  in  5°  steps. 
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IMPULSE:  A  FOkTkAN  CCCE  FOR  CALCULATING  THF  IMPULSE  RESPONSE  FOR 

A  bOTTOM-MOUNT £C  Ok  SUSPENDED  SOURCE  AT  A  SUSPENDED  RECEIVER.  FROM 
THE  RECIPROCITY  THtCkEM,  THIS  ALSO  SUFFICFS  FOR  CALCULATING  THE  IMPULS 
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this  routine  creates,  a  linear  f>iece"W1E£  continuous  function 

FITTING  The  DATA  IN  ARRAYS  t  ANC  X.  A  IS  REGARDED  AS  ThF 


Subroutine  —  PLTSSP 


SUBROUTINE  DMPcriCh) 


THIS  pCUTl\E  IS  *  PLUTlIfMt  A  £  J  U  N  C  T  TO  IMPULSE.  IT  USES  m  fiLF 
IN  WHICH  Afif  S  10k  E  l  t:  0  T  H  FLOTTJNG  LATA  AN  f1  INSTRUCTIONS.  1  Ht  STKUCTURt 


*E ACC10.10CG)A(1),T<I) 
CALL  LlNPLT (X , Y, TOC) 


READC10,1LCL)X(l).t(l) 
CALL  LINPlT (X ,Y , ’a  0) 
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APPENDIX  B.  LISTING  OF  TEST  PROGRAMS 


The  two  test  routines  presented  here  create  data  files  for  IMPULSE, 
first  creates  data  for  a  Hirsch  profile;  the  second  creates  data  for 
parabolic  mirror  bottom. 


Test:  THIS  totf  CHEATLS  A  DA*A  SFT  SL'I7APi.F  FOP  TEST  CASE  Fu 
Cf  IMPULSE,  AND  STOKES  THE  DATA  ON  A  FILE  (EXTFhnA(.  UNIT  >),  uSl 
ThF  HIRSCH  PRCFILl : 
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